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Abstract 

Properties of polyelectrolytes in tetravalent salt solutions are intensively investigated by a coarse- 
grained model. The concentration of salt and the size of tetravalent counterions are found playing 
a decisive role on chain properties. If the size of tetravalent counterions is compatible with the one 
of monomers, the chains show extended structures at low and at high salt concentrations, whereas 
at intermediate salt concentrations, they acquire compact and prolate structures. The swelling 
exponent of a chain against salt concentration behaves in an analogous way as the morphological 
quantities. Under certain condition, the electrostatics gives a negative contribution to the persis- 
tence length, in companion with a salt-induced mechanical instability of polyelectrolytes. Nearly at 
the same moment, it appears like-charge attraction between chains. The equal size of the tetrava- 
lent ions and the monomers is the optimal condition to attain the strongest attraction between 
chains and the most compact chain structure. Moreover, the ions form a multi-layer organization 
around a chain and, thus, the integrated charge distribution reveals an oscillatory behavior. The 
results suggest that charge inversion has no direct connection with redissolution of polyelectrolytes 
at high salt concentrations. 
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I. INTRODUCTION 



Solutions of highly-charged polyelectrolytes, such as biopolymers (DNA, proteins, etc.) 
and many synthetic polyelectrolytes, exhibit striking phenomena while multivalent salts or 
charged molecules are added into the solutions [1-4] . If the concentration of the added salt 
lies within a window of concentration, the solution is separated in two phases: one is a viscous 
phase containing condensed or precipitated polyelectrolytes, and the other is a dilute phase 
in which polyelectrolytes are dissolved. If it lies outside the window, the solution is in a 
homogeneous phase. These phenomena are called "reentrant condensation" [5] , or commonly 
referred as "DNA condensation" in the fields of biology and medicine [6]. This topic has 
attracted much attention recently due to its potential application in developing nonviral 
vectors of DNA for the purpose of gene therapy [7] . The amount of salt which results in the 
phase separation is determined by many factors, such as valence of salt, temperature, chain 
structure, and concentration of polyelectrolytes [2-4, 8]. It is known that the lower boundary 
of the window depends linearly on the concentration of polyelectrolytes but the upper one 
of the window is insensitive to it and roughly a constant. Both flexible polyelectrolytes and 
stiff polyelectrolytes display similar phase diagrams [3, 4, 9]. 

Despite decades of effort, our understanding on polyelectrolytes is still limited and many 
issues, including the reentrant condensation, are far from being resolved [10]. The aim of 
this paper is to investigate the properties of polyelectrolytes in salt solutions by means of 
computer simulations. There have been many simulations devoted to the study of poly- 
electrolytes [11-23]. However, most of them investigated the properties in salt- free solu- 
tions [11-16]. Only recently, due to the progress of computing power, large-scale simulations 
for systems with added salt can be realized. The works included, for example, flexible chains 
in monovalent or multivalent salt solutions [17-20], semiflexible chains with multivalent salt 
or charged molecules [21, 22], and rigid chains in the presence of multivalent salts [23]. But 
the salt concentration or the salt valence was not large enough to observe any significant 
evidence being able to connect to the redissolution of polyelectrolytes. Therefore, the reen- 
trant condensation could not be discussed thoroughly in these works. AUahyrarov et al. [24] 
recently studied the effective interaction between two immobile DNAs in the presence of 
tetravalent counterions. They have found that the effective interaction yields an attrac- 
tive minimum at short distances which then disappears in favor of a shallow minimum at 
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large separations, indicating the occurrence of DNA condensation and redissolution and a 
novel stable mesocrystal. In a more recent study [25], we have investigated the properties 
of mobile polyelectrolytes in tetravalent salt solutions and demonstrated that following a 
collapsed transition, polyelectrolytes undergo reexpanding transition upon addition of the 
salt. This phenomenon is just a single-chain version of the condensation and the redissolu- 
tion of polyelectrolytes. Moreover, by varying the ion size, we have shown that the excluded 
volume of ions plays a decisive role on the reentrant condensation. However, to well under- 
stand the phenomena, there is still a far way to go. For example, is it possible to obtain 
any information related to the phase boundary of reentrant condensation from simulations? 
There exist theories to explain the phase boundary [3-5] but it has never been studied in 
simulations. Concerning the size of polyelectrolyte in multivalent salts, Muthukumar and 
coworkers have applied single and double screening theories to explain the collapsed transi- 
tion of chains [18, 26, 27]. However, no existing theory can predict chain reexpansion in the 
region of high salt concentration. Moreover, chain morphology is essential to understand the 
properties of polyelectrolytes. It can be described by many quantities, such as asphericity 
and prolateness parameter [28]. To our knowledge, no such study has been reported yet for 
strongly-charged polyelectrolytes. An alternative way to understand polymers is to study 
scaling behavior. Although the concept of scaling has been successfully applied to neutral 
polymers [29], apphcation of this concept to polyelectrolytes [30-33] has encountered many 
difficulties owing to the long-range nature of Coulomb interaction. In this paper, we com- 
pute the single-chain structure factor. The swelling exponent is then calculated by studying 
the power-law-like regime of the structure factor. Since properties of polyelectrolytes depend 
strongly on ion size, we investigate in this paper the effect of ion size too. In our previous 
work [25], the sizes of ions were set to identical and varied simultaneously. As a consequence, 
Bjerrum association [34] was strongly intervened, particularly when ions were small. This 
setup has a disadvantage that the system was easily jammed by large ions, leading to the 
impossibility of investigation of the cases of big ions. Therefore, in this study we vary only 
the size of the multivalent counterions and fix unchanged the one of the other ions, so that 
the ion size effect is restrictedly investigated. 

In the study of polyelectrolytes, it is tempting to discuss the long range nature of Coulomb 
interaction using the concept of electrostatic persistence length. Odijk [35], Skolnick and 
Fixman [36] predicted that the electrostatic persistence length is proportional to the square 
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of the screening length rg. The deduction was based upon Debye-Hiickel theory. Conse- 
quently, the effective intramolecular interaction is always repulsive, resulting in a positive 
electrostatic persistence length [37, 38]. However, experiments have shown that electro- 
static persistence length can be negative under certain condition so that electrostatics gives 
a reversed contribution to chain flexibility [39, 40]. Simulations on this issue [41-45] were 
mainly performed within Debye-Hiickel approximation using a screened Coulomb poten- 
tial C/dh('^) oc e''^!'^" jr to describe the screening effect of salt without involving explicitly 
salt ions. A negative electrostatic persistence length, hence, has never been reproduced in 
the simulations. Noticeably, it is known that linear screening theory is valid only for a 
wcakly-chargcd system and breaks down for a highly-charged one, such as DNA and many 
biopolymcrs, due to strong ion correlations [17, 46, 47]. Recently, by employing a loop ex- 
pansion method which goes beyond mean field, Ariel and Andelman were able to predict the 
negative regime of the electrostatic persistence length for a strongly-charged rodlike poly- 
electrolyte [48]. Therefore, we expect that to observe a negative electrostatic persistence 
length, one should incorporate salt ions explicitly and use the bare Coulomb potential in 
simulations. Our study fulfills this requirement and hence offer a good chance to verify the 
above phenomenon and theory. 

Charge inversion is a universal phenomenon, occurred when a strongly charged macroion 
binds so many counterions that its net charge alters sign [49] . Nguyen ei al. intended to link 
charge inversion with the reentrant condensation of polyelectrolytes [5] . The idea was based 
upon that counterions form strongly correlated liquid on the surface of polyelectrolytes. 
In an intermediate salt region, the bare charges of the chains are almost neutralized by 
the condensed counterions; the correlation-induced attraction results in the condensation 
of polyelectrolytes. In a high-salt region, the counterions overcompensate the chain bare 
charges; Coulomb repulsion between chains dominates the attraction and, in consequence, 
the condensates are redissolved into the solution. Solis and Olvera de la Cruz gave a dif- 
ferent point of view using a theory called "two-state model" [50, 51]. They predicted that 
charge inversion is not necessarily occurred with redissolution of polyelectrolytes. Although 
ovcrcompensated by counterions, the redissolved polyelectrolytes can have either sign due to 
the association of colons [52]. This association depends strongly on the ion size which was 
not considered in the study of Nguyen et al. [5] . Therefore, Solis and Olvera de la Cruz pre- 
dicted a charge distribution which oscillates from a polyelectrolyte, in accordance with the 
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profile observed in our previous work [25] . In this paper, we investigate more extreme cases, 
including much smaller and much larger ions. We also study like-charge attraction between 
chains by calculating the potential of mean force and determine the conditions under which 
the attraction takes place. 

The rest of the article is organized as follows. Model and simulation details are given 
in Sec. II. Results are discussed in Sec. III. It is divided into five subsections. The first 
subsection (Sec. Ill A) is devoted to the study of the radius of gyration of a polyelectrolyte 
where salt concentration and polyelectrolyte concentration are varied. The results provide 
information helpful in understanding the boundaries of the condensation window. Sec. Ill B 
studies the morphology of multiple polyelectrolytes in salt solutions. We compute several 
quantities, including shape factors, asphericity, and prolateness parameter. The effect of size 
of multivalent ions is discussed. Single-chain structure factor is, then, calculated in Sec. Ill C. 
By least-square fitting the power-law-like regime, we compute the swelling exponent of a 
polyelectrolyte, which provides a fundamental vision how strongly-charged polyelectrolytes 
scale in salt solutions. Sec. Ill D deals with persistence length. The negative regime of the 
electrostatic persistence length is successfully reproduced in the simulations. Comparison 
of the results with theoretical predictions is made. Sec. IIIE provides a general picture of 
the effective interaction between polyelectrolytes by computing the potential of mean force. 
As long as knowing the condition to appear like-charge attraction between chains and the 
integrated charge distribution around a chain (Sec. IIIF), we are able to arbitrate if the 
phenomena of charge inversion and the redissolution of polyelectrolytes are related. We give 
our conclusions in Sec. IV. 



II. MODEL AND SIMULATION METHOD 



We employed a coarse-grained model to simulate a polyelectrolyte system. The system 
contains four bead-spring chains, one of which is composed of 48 beads (monomers). Each 
monomer dissociates a monovalent cation (counterion) into a solution and carries a negative 
unit charge. There are totally 192 monovalent cations in the solution. Solvent molecules 
are not modeled explicitly. Their effect is taken into account implicitly by considering them 
as a medium of constant dielectric constant. The geometry of the medium is cubic and 
periodic boundary condition is applied to simulate bulk environment. Salt molecules are 
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added into the system. In the solution, they are dissociated into cations (counterions) and 
anions (coions). We assume that the cations are tetravalent and the anions are monovalent. 
The amount of the cations and anions obeys the condition of charge neutrality. 

Three types of interactions are included in our simulations. The first one is the excluded 
volume interaction described by a purely repulsive Lennard- Jones (LJ) potential, 

, , / 4^Lj [{a/rr - {<T/rf] + £LJ for r < 
f^Lj(r) = < , (1) 

I for r > ^/2a 

where r is the distance between two particles, £lj is the interaction strength, and cr is the 
collision diameter. We assume that all the particles have identical £lj but the collision diam- 
eter may be different. The collision diameters for monomer, monovalent cation, tetravalent 
cation, and anion are denoted by dm, a+i, (T+4, and respectively. Lorentz-Berthelot 
mixing rule [53] is applied for the interaction between different kinds of particles. 

The second interaction is called finitely extensible nonlinear elastic (FENE) potential [54] 
and used to describe the bond connection between adjacent monomers on a chain. It reads 
as 

[^feneW = -2^FENE-Ro In (^1 - j^j , (2) 



where /cfene is the spring constant, ^ is bond length, and Rq is the maximum extension of a 
bond. All particles interact via the third interaction. Coulomb interaction, which is written 
as 

C/coui(r) = A;bT4^ (3) 

where ZjC is the charge of the ith particle, T is the temperature, k-Q is the Boltzmann 
constant, and £b is the Bjerrum length defined as the separation distance at which the 
Coulomb potential between two unit charges e equals to the thermal energy k-oT. £b is 
equal to e^/ (ATreeokBT) where e is the dielectric constant and eo is the vacuum permittivity. 

We perform molecular dynamics (MD) simulations in this study [75]. We set cr+i and 
equal to and vary the size of tetravalent cation o"+4 from O.Oa^ to 4.0(Tm to investigate 
the effect of ion size on properties of the polyelectrolyte solution. Simulation parameters are 
chosen as follows: £lj = 0.8333/i;Br, /cfene = 5.8333A;b2^/c^, and Ro — lo^. The Bjerrum 
length is set to ^b = ScTm- Similar systems in salt- free solutions have been investigated by 
Stevens and Kremer [14]. We assume that the masses of all the particles are identical and 
equal to m, and choose and k-^T to be the length unit and the energy unit, respectively. 
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Therefore, the natural time unit is t — 



Cm^/^fn/Jk^T). The system is studied in canonical 



ensemble where the temperature is controlled using Nose-Hoover thermostat with the fre- 
quency of heat transfer with thermostat equal to l.Or"^. Coulomb interaction is calculated 
by the technique of Ewald sum. Except in Sec. Ill A, the monomer concentration is fixed at 
Cm = 0.008(7"^ which corresponds to a box size L = 28.8(Tm, whereas the salt concentration 
is varied. The amount of salt is controlled such that the volume fraction of total particles 
is smaller than 30%. Equations of motion are integrated applying Verlet algorithm where 
time step is set to O.OOSr. Each simulation starts with an equihbration phase of 2 x 10^ MD 
steps and is followed by a production run of 10^ MD steps in which data are collected every 
1000 steps. In order to shorten the notation, we will use Cm as the length unit, cx"^ as the 
concentration unit , and e as the charge unit in the following text. 

III. RESULTS AND DISCUSSIONS 

A. Radius of gyration of a polyelectrolyte at different monomer concentrations 

We started our study by simulating a single polyelectrolyte with added tetravalent salt 
at various monomer concentrations. The chain consists of 48 monomers and the size of the 
tetravalent cations is set to be identical to the one of the monomers. We used the radius 
of gyration to characterize the chain size. The radius of gyration is calculated by the 
equation 



where N is the number of monomers on the chain, is the position vector of the ith 
monomer, and fcm the center of mass of the chain. The results are presented in Fig. 1 where 
the mean square radius of gyration (R^) is plotted as a function of salt concentration Cg at 
four monomer concentrations. Cm = 0.00032, 0.002, 0.008 and 0.016. The angular bracket 
(■) means averaging over chain configurations. Two behaviors were observed. The first one 
is that each curve displays a similar trend of behavior against salt concentration. Upon 
addition of salt, the chain undergoes two structural transitions: a collapsed transition and 
a followed reexpanding transition, as having been observed in the previous work [25]. The 
first transition results from the screening of the Coulomb interaction by the added salt. The 
second transition is still not well understood although it is a result of the competition between 
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FIG. 1: (i?g) as a function of Cg at four monomer concentrations, Cm = 0.00032, 0.002, 0.008, 
and 0.016. The symbols used to represent the four Cm are indicated at the left-bottom side of the 
figure. The value of {Rg) in a salt-free solution (Cs=0.0) is depicted near the left vertical axis. The 
four curves in the regime Cs < C* can overlap each other (shown in the inset) under a suitable 
transformation described in the text. 



energetic gain and entropic loose in counterion condensation and in shape transformation of 
polyelectrolytes to minimize the free energy. We observed that the minimum of {R^) appears 
near the salt concentration C* = Cm/4 and the more dilute the monomer concentration, the 
smaller value the minimum of (i?g) will be. C* is the equivalence point at which the amount 
of charges of the tetravalent counterions exactly neutralizes the charges on polyelectrolytes 
as if the salts and the polyelectrolytes are in the stoichiometric reaction. We will call C* 
the equivalence concentration. In a salt-free solution {Cg = 0), we found that the chain 
size increases with decreasing monomer concentration (refer to the four data points near 
the left vertical axis of the figure). This effect has been discussed in literatures [14, 33] and 
our results are in accordance with them. And the second behavior is that (R^) for different 
Cm attains approximately the same value while the salt concentration is larger than the 
equivalence concentration. Therefore, the curves approximately overlap each other. We 
noticed, furthermore, that in the region Cg < C*, the (i?g) curves for different Cm can also 
lie over each other under an appropriate transformation. We demonstrate this overlap in 
the inset of Fig. 1 where the abscissa is Cs/Cm and the ordinate is A^{{R^) — (-R*^)) with 
(i?*^) the minimum value of (-Rg). Am depends smoothly on Cm and, in the figure, is equal 
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to 1, 1.267, 1.657, and 2.049 for = 0.00032, 0.002, 0.008, and 0.016, respectively. 

The size of polyelectrolytes in multivalent salt solutions has been studied by small-angle 
scattering experiments [55]. The results showed that it diminishes with increasing or de- 
creasing Cs toward the region where the condensation of polyelectrolytes occurs (cf. Table 
2 of Ref. [55]). The results of our simulations are in accordance with the experiments. We 
make a hypothesis that chain size can be served as a qualitative condition to distinguish 
between a monophasic solution of polyelectrolyte and a diphasic one. More precisely, if the 
chain size is larger than some value, the chains are dissolved in the solution and the system 
is in a homogeneous phase; if it is smaller, the chains are precipitated from the solution and 
a phase separation occurs. Under this hypothesis, polyelectrolyte solutions of different 
become monophasic as Cs is increased over a concentration in the high-salt region, since the 
chain size exceeds the some value nearly at the same Cg due to the similarity of the (R^) 
curves against Cg in this region. Therefore, the upper boundary of the salt concentration 
window for the condensation of polyelectrolyte is independent of Cm. Similarly, the lower 
boundary of the window approximately depends linearly on Cm, owing to the similarity of 
the ylm((-Rg) — (-Rg^)) curves against Cs/Cm in the low-salt region; the weak dependence of 
Am and (R*^) on Cm ensures that the chain size surpasses the some value roughly at the same 
Cs/Cm. Therefore, Fig. 1 provides an essential information to understand the boundaries of 
reentrant condensation. 

The aim of this article is to understand the properties of multiple flexible polyelectrolytes 
with added salt in an aqueous solution. It is known that the Bjerrum length in water is 
equal to 7.14A at room temperature. Our simulation setup, hence, corresponds to a collision 
diameter dm of monomer equal to 2.38A. This model represents a polyelectrolyte system 
with linear charge density equal to a unit charge per 2.62A[76], and can be employed to 
simulate a strongly-charged system such as sodium polystyrene sulfonate. The four monomer 
concentrations studied in this section, while converted to real unit, are equal to 0.039M, 
0.246M, 0.986M, and 1.971M, respectively. Even the largest Cm is still smaller than the 
overlap threshold estimated by 3A^/(47ri?g). Therefore, our systems are dilute polymer 
solutions. Readers should be aware that the value up to nearly 2M is a very high monomer 
concentration for a real polymer system, and our systems are dilute only in the formal 
meaning of the world since the chains are very short. In the following sections, we fix Cm 
at 0.008 and simulate a system comprising multiple chains. This monomer concentration 
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is one or two order of magnitude higher than a typical monomer concentration used in 
experiments [3, 4, 55]. However, based on the fact that systems at different monomer 
concentrations reveal similar evolution of chain size against salt concentration, the results 
obtained at Cm = 0.008 can be used to understand a system at more dilute Cm. We mention 
that this Cm is a choice so that the simulation box is not so big. Therefore, the amount of 
the added salt is controllable under limited computing resources. Consequently, to study the 
behavior of polyelectrolytes in salt solutions, covering a broad range of salt concentration, 
becomes numerically feasible. 

B. Shape of polyelectrolytes 

In this section and in the followings, the studied system contains four polyelectrolytes 
and the monomer concentration is fixed at Cm = 0.008. This section is devoted to the study 
of chain conformation in tetravalent salt solutions. The effects of salt concentration and size 
of tetravalent counterions are investigated. We vary the size of tetravalent counterion, (7+4, 
from 0.0 to 4.0 and fix the size of the other ions at 1.0. A large a +4 can represent a large 
cation, a large ion group, or a large charged colloid. It can be also a result of the hydration 
of a small ion. A small (T+4 may be less representative for real systems but theorists are 
generally interested in it because many models are established on the base of vanishing ion 
size. 

Morphology of a polymer can be quantified by many quantities. The first quantity that 
we studied is the shape factor rj, defined as the ratio of the mean-square end-to-end distance, 
(Rl) = ((ri — rAf)^), to the mean-square radius of gyration (-Rg). r] attains a value 12 for 
a rodlike polymer, 6.3 for a flexible chain in a good solvent, and 6 for an ideal chain [14]. 
For a compact structure, the value of rj is small. For example, if the chain has a spherical 
morphology and the two ends are randomly positioned, r] is equal to 2. In Fig. 2 (a), we 
show how the shape factor 77 varies with Cg at different (j+4. In the figure, a reference 
curve, marked by "neutral", is depicted. It represents a system composed of four neutral 
polymers of 48 monomers in tetravalent salt ((j+4 = 1.0) solutions. Focus firstly on the 
polyelectrolyte solutions with (T+4 = 1.0, 1.4 and 2.0. We witnessed that r] decreases at first 
and then increases upon addition of salt, in analogy of the radius of gyration studied in the 
previous section, and shows a V-shaped curve against Cg in the semilog plot. In the low-salt 
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FIG. 2: (a) r? = {Rl)/{Rl) as a function of Cs, and (b) ^ = {Riy/^/{Rh) as a function of C^, 
for the cases of different size cr+4 of the tetravalent counterions. The symbol used for each (T+4 
is indicated at the left-bottom side of the figures. The error bars are smaller than the size of the 
symbols. The result for neutral polymers is depicted as reference and marked by "neutral". 



region or in the high-salt region, rj is larger than or close to the value 6.3, indicating that 
the chains have elongated or coil-like structures. In between, it appears a minimum near 
the equivalence concentration C* = 0.002. The value of 77 tells us that the chains possess 
compact structures but not so compact as spheres. For the cases with larger (T+4 (cr+4 = 3.0 
and 4.0), 1] shows a shallow minimum and the whole curve stays basically above the value 
6.3. In these cases, the chains are in an extended state. It is worth to notice that for 
(T+4 > 1.0, T], in principle, increases with (T+4 at a given Cs- Now pay attention to the cases 
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with (7+4 < 1.0. 77 displays a similar V-shaped curve against Cg but opposite to the previous 
cases, its value increases with decreasing (7+4. Moreover, the degree of chain reexpansion in 
the high-salt region is weakened while 17+4 is decreased. Noticeably, for the case 17+4 = 0.0, 
1] displays a sharp transition from a value 7.2 in the low-salt solution to a value 6.5 in the 
vicinity of C*. In the latter region, the chains behave similarly to neutral polymers. 

The second quantity that we studied is ^, which is the ratio of the radius of gyration 
(i?2)V2 to the hydrodynamic radius (R^) computed by R^' = N-^^ti Ef=ij^i^^')- C is 
equal to ~ 0.775 for a hard sphere and attains a value 2.25 for a rodhke chain [56]. 

Renormahzation group theory has predicted ^ ~ 1.56 in a good solvent and ^ ~ 1.24 at © 
point [57], both of which have been shown in agreement with experiments [58]. We present 
in Fig. 2(b) the variation of ^ against Cg- The referenced curve for the system of neutral 
polymers is depicted. The curve was found very close to the value 1.24, which seems to 
indicate an ideal-chain behavior. However, the neutral chains should behave as random 
coils and the value of ^ would be larger than 1.24. The deviation of ^ from the value of a 
coil is probably due to the short chain length used in our simulations. Diinweg et al. have 
studied the finite size effect and shown that the dependence of ^ on chain length is not 
negligible [59]. Therefore, ^ reported here departs slightly from the asymptotic value of an 
infinite chain length. Nonetheless, the results are still indicative and useful to understand 
how ^ varies with salt concentration. In the small Cg region, the polyelectrolytes are more 
extended than the neutral polymers. In the middle region of Cg, the chains attain compact 
structures when 0.5 < (7+4 < 2.0. The dependence of ^ upon (7+4 resembles to that of rj. 
For large cr+4 ((7+4 > 2.0), the chains are elongated. For vanishing (7+4, we observed that 
following a sharp decrease, ^ reattains approximately a constant value, close to that for the 
neutral polymers. 

A more fundamental way to study the shape of a polymer is to investigate the tensor of 
radius of gyration T [60] , defined by: 

^ N N 

= ^ ~ ^J") (^^^3 ~ '^^■Z^) 

i=l j=l 

where a, (3 — 1,2,3 denote the three Cartesian components and ria is the a component of the 
position of the ith monomer. Let Ai, A2 and A3 be the three eigenvalues of T. A quantity 
called "asphericity" [28], which measures the deformation from a spherical geometry, is 
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defined by 

1 / (Ai - \2f + (A2 - + (A3 - Ai)^ 

2 \ (Ai + A2 + A3)2 
It takes a value between (sphere) and 1 (rod). For a coil polymer, it is 0.431 obtained 

from simulations [61]. Another quantity describing the degree of prolateness is defined by 

/ (Ai-A)(A2-A)(A3-A) 
\ A3 

where A = (Ai + A2 + A3)/3 [28]. It is ranged from —1/4 to 2 and can be used to distinguish 
between a prolate shape (0 < P < 2) and an oblate one (— | < P < 0). For neutral polymers 
in the dilute limit, simulations predict a prolate chain with P = 0.541 [62] . We present in 
Figs. 3(a) and 3(b) the asphericity A and the prolateness parameter P as a function of Cg. 
Resembling to the behaviors of r] and ^, A and P show a minimum value in the vicinity of C*. 
The curve associated to (7+4 = 1.0 roughly separates two different behaviors of evolution: 
A and P increase with increasing and with decreasing (7+4 from 1.0. The value of A lies 
between 0.275 and 0.55. Therefore, even in a compact structure (where Cg is near Cg*), the 
chain shape is asymmetric. In a low-salt solution, the chains show extended morphology 
but are not as elongated as rods. On the other hand, P is larger than zero for all of the 
studied cases, indicating that the polyelectrolytes display prolate conformation. At a high 
salt concentration, A is close to 0.45 and P is close to 0.54, suggesting once more that the 
polyelectrolytes behave as the neutral polymers in the high-salt region. For the case with 
null (T+4, chain reexpansion does not occur; A and P stay at a constant value when Cg 
goes beyond C*. The results obtained in this section strongly suggest that 17+4 = 1.0 is the 
optimal condition to pack polyelectrolytes into the smallest volume. 





C. Single-chain structure factor and swelling exponent 

In this section, we investigate the single-chain structure factor. This quantity is ex- 
perimentally accessible and provides an essential information to theoretical analysis. The 
single-chain structure factor is defined by 



1 



iVp 



N 



^exp(ig-fpj) 



(8) 



where q is the scattering wavevector, A^p = 4 is the number of polymer chains. A" = 48 is 
the number of monomers on a chain, and r^j is the position of the jth monomer on the pth 
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FIG. 3: (a) A as a function of Cg, and (b) P as a function of Cg, for different cr+4. The symbol 
used for each (T+4 is indicated at the left-bottom side of the figures. The error bars are smaller 
than the size of the symbols. The result for neutral polymers is depicted as reference and marked 
by "neutral". 



chain. Suppose that the system has spherical symmetry. The dependence of the structure 
factor on the orientation of q can be therefore integrated out: S{q) = ^ J dQS{q) where 
q is the norm of q and Q is the solid angle. This integration was performed numerically 
by sampling over a set of randomly-oriented q vectors. We present, in Fig. 4(a), S{q) at 
several salt concentrations for the case of a+4 = 1.0. For the reason of clarity, the curves 
in the figure have been multiplied by 2" with n being increased by one for each curve from 
the bottom to the top. S{q) for the neutral polymers is depicted as reference with n = 0. 
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We witnessed that S{q) shows typical behavior. In the Guinier regime {qRg <^ 1), S{q) 
behaves similarly to A'"(l — q^{R^)/3). It is demonstrated by plotting A'"(l — q'^{R^)/3) (in 
dotted curve) near the corresponding S{q) with the value of (i?g) adopted from Sec. IIIB. 
The dotted curve terminates at g = R^^ in order to indicate the boundary of the Guinier 
regime. In the regime R^^ <^ g -C cr~^, S{q) shows power- law-like behavior, manifested by 
a linear dependence of logS'(g) on logg. The slope s of the linear dependence is related to 
the swelling exponent u by the relation 1/ — —1/s, which describes the size extension of a 
polymer: R^ ~ N'^. In the regime q > cr"^, S{q) displays non-universal behavior. However, 
while q ^ (7^^^, self-scattering of monomer is the only contribution, leading, therefore, S{q) 
to a universal value 1. We mention that the Krathy regime, in which S{q) ~ q^^, is not 
observed in our study because the system is a dilute solution and the chain length is not 
long enough [63]. For (7+4 other than 1.0, we observed analogous trends of behavior. For 
example, the variation of S{q) for different (7+4 at the equivalence concentration is presented 
in Fig. 4(b), in which we rediscovered the three regimes: the Guinier regime, the power- 
law-like regime and the non-universal regime. We noticed that in the last regime, a small 
bounce appears on the curve and, with increasing the size of tetravalent ions, gradually 
moves toward small q. 

The swelling exponent u of the polyelectrolytes was calculated by performing least-square 
fits in the log-log plot of S{q) in the power-law-like regime. The fitting region was chosen 
to be 0.63 < q < 1 and the fitting rehabihty has been verified to be excellent. We present 
the results of the swelling exponent i/ in Fig. 5. 1/ for the referenced neutral polymers was 
verified at first and found to attain a value close to the Flory's one, 0.6. Next, we focused on 
the swelling exponent for the polyelectrolytes. We observed that u exceeds the value for the 
neutral polymers at a low salt concentration. It indicates a strong extension of chain size due 
to Goulomb repulsion between monomers. For the cases with (7+4 = 0.5, 1,0, 1,4, and 2.0, 
the swelling exponent shows a V-shaped curve against Cg: u decreases and then increases 
with increasing salt concentration. A minimum occurs near the equivalence concentration 
C* and, at this moment, the chains show compact structures but are less compact than 
spheres because a spherical structure scales as N^^^. We remark that the swelling exponent 
obtained here for (T_|_4 = 1.0 is consistent with that obtained by direct variation of the chain 
length of a system containing a single polyelectrolyte [64] except in the mid-salt region where 
u in the latter study is slightly smaller and acquires a value 1/3 at C*. This deviation is 
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FIG. 4: (a) S{q) for a+4 = 1-0 at different Cs and (b) S{q) at = 0.002 for different cT+4. The 
value of Cs or cr+4, at which the simulations were running, is indicated at the left side of each 
corresponding curve. S{q) for the neutral polymers at Cs = 0.002 is depicted as reference and 
marked by "neutral". For clarity, the curves, from bottom to top, have been multiplied by a factor 
2" with n started from zero and increased by one for every curve. The dotted curve, near each S{q) 
and terminated at q = Rg^, is the Guinier function A^(l — q^(i?g)/3). Two dashed lines, which 
represent two power-law dependences, q^^ and q^'^ (marked by and "-2", respectively), are 
drawn in the plots. 
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FIG. 5: as a function of Cg for different cr+4. The symbol used for each fT+4 is indicated at the 
left-bottom side of the figure, for the referenced neutral-polymer system is depicted and marked 
by "neutral". 

due to multi-chain aggregation occurred in our system, which deforms the chain shape from 
a sphere and, consequently, increases u. For the case with vanishing (T+4, following a drastic 
decrease, u retains roughly a value 0.6 in the high-salt region and the chains, hence, scale 
like the neutral polymers. For the cases with large cr+4 (cr+4 = 3.0 and 4.0), u is greater 
than 0.6 but still less than 1.0. The chains, hence, do not reveal rigid rod-like structures. 
In principle, (T+4 = 1.0 separates two behaviors of evolution at a salt concentration: u 
increases with increasing and with decreasing (T+4 from 1.0. The results demonstrate again 
that the polyelectrolytes are mostly condensed while the size of the tetravalent counterions 
is compatible with the one of the monomers. 

D. Persistence length 

There is a considerable interest to understand the variation of the persistence length of a 
polyelectrolyte against salt concentration. Persistence length is a local property, which mea- 
sures the stiffness of a polymer, and is distinguished from global conformational quantities. 
There are many ways to define the persistence length [42, 65]. In this section, we adopt a 
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microscopic definition to calculate it: 

^ (iV/2)-l 



2b ^ 

i=0 



{bN/2 ■ b(N/2)-i + bN/2 ■ b(^N/2)+i) 



(9) 



where bj = rj+i — Vj is the jth bond vector and b = (&|)^'^^ is the root-mean-square bond 
length. In this definition, the contribution of the finite bond length, which corresponds to 
the 2 = term in Eq. 9, is considered. The obtaining £p as a function of Cg for different 
cr+4 is plotted in Fig. 6. Resembling to the quantities of chain conformation, £p decreases 
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FIG. 6: Cp as a function of Cg for different cT+4. The symbol used for each cr+4 is indicated at the 
left-bottom side of the figure. The error bars are smaller than the size of the symbols. The bare 
persistence length £0 is depicted as reference and marked by "neutral" . 



in the chain collapse region and increases in the chain reexpansion region. A minimum 
value appears in the vicinity of C*. The degree of increase diminishes with decreasing cT+4. 
For vanishing cr+4, it is so weak that £p is roughly a constant as > C*. Similar to the 
results in the previous sections, ip increases with both increasing and decreasing a+4 from 
1.0. cr+4 = 1.0, therefore, denotes a condition to separate the two behaviors of evolution of 

The persistence length has generally two contributions: the bare persistence length (in 
the absence of electrostatic force) and the electrostatic persistence length (owing to the 
Coulomb interaction) [35, 36]. In this study, the bare persistence length was calculated by 
performing simulations on the referenced system where four neutral chains and 192 neutral 
particles are placed in salt solutions. The result is depicted in Fig. 6 and the curve is marked 
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by "neutral" . We found that £o is not sensitive to the salt concentration and roughly attains 
a value 2.3. It shghtly decreases in the high-salt region, probably due to the jamming effect 
by the salt ions. For the polyelectrolyte solutions with large (T+4 (cr+4 = 3.0 and 4.0), the ip 
curve thoroughly lies above the £0 curve, indicating a more rigid flexibility than a neutral 
polymer. The electrostatic persistence length is, hence, positive. On the other hand, for 
the polyelectrolyte solutions with (7+4 < 2.0, ip is smaller than £q in the mid-salt region and 
shows one or two intersections with the io curve. Therefore, ie — ip — io is positive in a 
salt-free solution and becomes eventually negative upon addition of salt. It increases in the 
high-salt region and reattains a positive value for 17+4 = 1.4 or 2.0. 

The negative contribution of electrostatics to the persistence length has been experimen- 
tally observed [39, 40]. Recently Ariel and Andelman (AA) pointed out that a negative 
indicates a mechanical instability (collapse) of polyelectrolytes due to the presence of mul- 
tivalent counterions [48]. By taking into account the thermal fluctuations and correlations 
between bound counterions, they derived an expression of ie for an intrinsically-rigid charged 
polymer in (Z:l)-salt solutions in the hmit of 6 <S k"^ <^ {N — 1)6: 

where F = Zis/b is the Coulomb strength parameter between a monomer and a Z-valent 
counterion, k — y/4:7rZ{Z + 1)£bCs is the inverse Debye-Hiickel screening length, and i^^^ — 
{AZ"^ i^)"^ is the electrostatic persistence length deduced by Odijk-Skolnick-Fixman (OSF) 
theory [35, 36]. In the derivation, they have assumed point-like charges for the ions. We 
present, in Fig. 7, i^ obtained from the simulations and that predicted by Eq. (10) as a 
function of k. We witnessed that i^ depends strongly on the size of the tetravalent ions. 
The negative regime of i^ has been successfully captured in the simulations. To make a 
comparison with the theoretical prediction, we focused on the case with vanishing (7+4. 
Please notice that this case does not exactly fulfill the assumption of the theory because 
in the simulations, the monovalent counterions and colons are not point charge-like and, 
also, the chains are not rigid. However, according to our experiences [25], ions other than 
multivalent counterions do not play a decisive role on the properties of polyelectrolytes. 
This comparison, therefore, makes sense if the effect of chain rigidness is negligible. We 
emphasize that the electrostatic correlations are fully taken into account in the simulations. 
In our study, the region in which Eq. (10) is vahd is 0.02 <^ k <^ 0.91 because Z — A, i-Q — 
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FIG. 7: ic as a function of k for different cT+4. The symbol used for each cj+4 is indicated at the 
left-bottom side of the figure. The error bars are smaller than the size of the symbols. The dashed 
curve depicts ie predicted by Eq. (10). 

and b is about 1.1. Notice that Eq. (10) has a singularity at k = 1/6 ~ 0.91 and, hence, 
the curve is divided into two branches. The left branch is located in the valid region but 
fails to predict ie- It gives not only a wrong value but also a wrong trend of variation of ie 
against k. We noticed that Eq. (10) is a modification of i^^^ multiplied by a factor but the 
valid region given by AA theory [48] is incompatible with that of OSF theory. OSF theory 
is held when the screening length is much smaller than the bare persistence length [44]. 
Therefore, AA theory may be true in the region n ^ 1/b. As shown in the figure, ie for 
vanishing (T+4 is well described by the right branch of the curve of Eq. (10), which seems to 
support this conjecture. This topic deserves a further analysis in the future. 

E. Potential of mean force 

Experiments and simulations have shown that like-charged macroions can attract each 
other in a solution [66-68]. It is hence relevant to know under which conditions the like- 
charge attraction takes place between polyelectrolytes. To answer this question, we investi- 
gate "potential of mean force" in this section. Potential of mean force describes the effective 
interaction between two particles (molecules) in a medium, which includes two contribu- 
tions: the direct interaction between two particles (molecules) and the indirect interaction 
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via other particles (molecules) and the medium. This quantity can be calculated by simu- 
lations: Wapir) — —kBTln Qa/sir) where gap{f) is the radial distribution function between 
particles of species a and gapi^f) is computed by 

9 sir) = (11) 

where and A?^ are the number of a-particles and that of /3-particles, respectively, V is the 
volume of the simulation box, and Ha/sir) — '^f^i hipir) with hipir) denoting the number 
of /9-particles lying in a spherical shell, centered at the ct-particle i, of radius r and shell 

thickness 5r. 

Because our system is composed of only four chains, it is not easy to reduce statistical error 
if we directly calculate the potential of mean force between chains. We, therefore, calculated 
the potential of mean force between monomers on different chains, which also provides the 
information about the chain interaction. We denote this quantity VFmm'('") ^^'^ the results 
are presented in Fig. 8 which contains six plots, from top to bottom and from left to right, 
corresponding, in orders, to the cases (j+4 = 0.0, 0.5, 1.0, 2.0, 3.0 and 4.0. PFmmK^) ^^e 
referenced neutral polymers is also depicted in the figure. In a salt-free solution, H^mm'(^) 
displays a repulsive interaction, and is stronger than that for the neutral polymers owing to 
the electrostatic repulsion between monomers. This can be seen by comparing the range of 
repulsion estimated by equaling W- m-m 'ir) to the thermal energy k-oT. The range of repulsion 
is at r = 7.1 for the polyelectrolytes and at r = 3.0 for the neutral polymers. Upon addition 
of salt, VFjiimK^) responses in two distinct ways, which depends on the ion size. The first one 
is applied for small or large tetravalent counterions, such as cr+4 = 0.0, 0.5, or 4.0. In these 
cases, Wm.m'i'r) is purely repulsive in the course of addition of salt. Hence, the monomers on 
different chains repel between themselves and the system exhibits no like-charge attraction. 
The second way is happened for an intermediate size such as (7+4 = 1.0, 2.0, or 3.0. At 
this moment, VFmmK^) displays an attractive well in a mid-salt region. Therefore, it appears 
like-charge attraction between chains. In the first case, the range of repulsion of Wmin'(^) 
decreases to a minimum value nearly at Cg = 0.002. Further increase of Cg either leaves 
M^mm'(^) unchanged (for vanishing (T+4 and for (T+4 = 4.0) or increases the range of repulsion 
of Wrara'{r) (for (7+4 = 0.5). Wc wiU call Wra-m>{r) at Cs — 0.002 "the curve at equivalence", 
which separates two behaviors of evolution of VFmm' {f) as Cg is increased or decreased from 
Cg*. We noticed that Wy^^f{r) tends toward the referenced curve of the neutral polymers at 
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FIG. 8: VFmm'(r) for cj+4 = 0.0, 0.5, 1.0, 2.0, 3.0 and 4.0, shown, in orders, from top to bottom 
and from left to right, in six plots. Each plot includes several Wmm/(r) at different Cg. The value 
of Cs is indicated in the plot or near the corresponding curve. Wnim'(^) neutral polymer 

system is depicted as reference and marked by "neutral". The unit of ordinate axis is /cbT. 



high sah concentrations except for the deviation happened at small r due to the excluded 
volume of large ions. For vanishing cr+4, the curve at equivalence is very close to the reference 
curve and Wuim'{r) is roughly unchanged while Cs > 0.002. For (T+4 = 0.5, the curve at 
equivalence lies below the referenced curve; Wmm'{r) is bounded between the two curves 
and increases with Cs- And in the second case, VFmm'(^) displays a potential well located 
near r = 0"^ + a ^4. It suggests an attraction mediated by a tetravalent counterion. The 
depth of the well can be more profound than fceT. Therefore, a stable bound state can be 
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formed between monomers on different chains and it is a prerequisite to occur salt-induced 
phase separation. The maximum depth of the well is happened while Cg is C*. In general, 
we found that W^-ni'{r) is bounded between the salt-free curve and the curve at equivalence 
while Cs < 0.002, and between the latter curve and the referenced neutral curve while 
Cs > 0.002. 

We have shown that the occurrence of the like-charge attraction depends on both salt 
concentration and ion size. In Fig. 9 we present, in the variable space, the region where 
W rr,m '(r) shows an attractive well. In the figure, the data points (symbolized by star) denote 




c 

s 

FIG. 9: Two distinct regions in the (Cg, cr-|_4) -variable space. In Region (I), W^^'{r) shows an 
attractive well and in Region (II), it is purely repulsive. 

the salt concentrations for a given cr+4, at which a purely repulsive Wmm'('") turns to show 
an attractive well and vice versa, and the closed solid curve is sketched to help readers to 
separate the two regions. It clearly shows that an effective attraction appears only when (j_|_4 
is neither too small nor too big and salt concentration is intermediate around the equivalence 
concentration. We mention that the occurrence of an attractive region in the potential of 
mean force is not sufficient to have a phase separation of the system. Phase transitions 
in charged colloid systems have been largely discussed in literatures but there have been 
relatively few studies of phase transitions in charged chain systems (see Refs. [69] and [70] 
and references therein). 

Fig. 10 shows snapshots of the simulations for (7+4 — 0.5, 1.0 and 4.0 at two salt concen- 
trations Cs — 0.002 and Cg = 0.008. At Cg = 0.002, we saw that the polyelectrolytes favor 
to be apart from each other for both of the cases (7+4 = 0.5 and (7+4 = 4.0. However, delicate 
difference was witnessed that the previous case displays a single-chain collapse but the latter 
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(c) 

FIG. 10: (Color online) Snapshots of simulations for (a) (J4.4 = 0.5, (b) cr+4 = 1.0, and (c) 
cr+4 = 4.0, at Cg = 0.002 (left) and Cg = 0.008 (right). In the pictures, the bead-spring chains are 
represented in gray color (red color online) and the monovalent counterions in light one (white color 
online). The tetravalent counterions are represented in dark-colored (blue-colored online) particles 
whereas the colons are not plotted for the clarity of the pictures. In the snapshots, a chain may 
appear more than one fragment because of pei^^dic boundary condition. 



one does not. For a+4 — 1.0, the chains favor to aggregate together, in accordance with the 
results of Fig. 8 which show the occurrence of an attraction between the negatively-charged 
polymers. This attraction is mediated by counterions, mainly the tetravalent counterions, 
since almost all the tetravalent counterions are condensed on the polyelectrolytes. We point 
out that the condensation of tetravalent counterion is not a sufficient condition to occur 
chain aggregation. It can be seen from the case (7+4 = 0.5 in which most of the tetravalent 
counterions are condensed but the chains does not aggregate. At the high salt concentra- 
tion Cg = 0.008, the attraction between chains is reduced, in comparison with the cases at 
Cg 0.002. The chains segregate or show tendency of segregation from each other. In the 
simulations, we observed that while Cg > 0.002, only part of the tetravalent counterions can 
condense onto the chains. 

F. Integrated charge distribution around a chain 

In this section, we investigate the charge distribution around a polyelectrolyte. It en- 
ables to understand how ions distribute around a chain and provides useful information for 
theorists to develop models. We define a worm-like tube of radius r around a chain to be 
a union of spheres of radius r centered at each monomer on the chain. We calculated 
the averaged total charge inside a worm-like tube, including the charge of the ions and the 
one of the monomers. This quantity is called "the integrated charge distribution" and is 
denoted by Q{r). 

We show in Fig. 11 Q{r) at various Cg for (7+4 ranged between 0.0 and 4.0. Focus firstly 
on the case with vanishing (7+4. In the region r < 0.5, Q{r) attains a value equal to the 
bare charge of a chain, —48. A stair-like behavior then appears for Q{r) in the region 
0.5 < r < 1.0 while 0.0 < Cg < 0.002. Outside the region, Q{r) monotonically increases 
and tends toward zero at large r due to electroneutrality. Since r < 1.0 is the depletion 
zone for the particles other than tetravalent counterions, the behavior of Q{r) in this region 
is completely determined by the condensation of tetravalent counterions. The stair-like 
behavior indicates a perfect condensation on the surface of the chain, owing to the vanishing 
size of the tetravalent counterions. For instance, at Cg = 0.001, we have 24 tetravalent 
counterions in the simulation box and, therefore, 6 condensed tetravalent counterions on 
a chain, in average. A perfect condensation will lead to a value —24 for Q{r) and this is 
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FIG. 11: Q{r) at different Cg for c7+4=0.0, 0.5, 1.0, 2.0, 3.0 and 4.0, shown, in orders, from top to 
bottom and from left to right, in six figm'es. The Hne arts used for different Cg are indicated at 
the right-bottom side of each figure. 



exactly what we have observed in the figure. At the equivalence concentration, Q{r) is zero 
while r > 0.5. If Cs is further increased (Cg > 0.002), Q{r) reveals a positive peak in the 
region 0.5 < r < 1.0, which indicates an overcharging, and, then, rapidly decays to zero as 
r > 1.0. The hight of the peak suggests that only a portion of the tetravalent counterions 
are condensed on the chain. For the case (T+4 = 0.5, the region in which only the tetravalent 
counterions can intervene is narrowed to 0.75 < r < 1.0. In the region, the step of the 
stair-like Q{r) becomes tilted while Cg < 0.002. If Cs > 0.002, Q{r) displays a positive 



26 



peak, which is more pronounced than that for vanishing (7+4. For the case a +4 — 1.0, an 
oscillatory behavior occurs at high Cg. Similar phenomena have been observed in other 
systems such as charged plate [71], colloid system [72], and rigid polyelectrolyte [23]. This 
oscillatory behavior indicates a multi-layer organization of ions: in the close neighborhood 
of a polyelectrolyte, tetravalent counterions condense and form a positively-charged layer 
which overcharges the chain; exterior to this layer, anionic particles intervene and form a 
negatively-charged layer which overcompensates the charge inside and, hence, the integrated 
charge becomes negative, and so forth. Notice that the depletion zone around a chain for 
tetravalent counterions is r < (1 -|- 17+4 )/2. It becomes larger than that for the other ions 
while a_|_4 > 1.0. In these cases, Q{r) decreases with increasing Cg in the depletion zone, but 
shows completely opposed behavior just outside the zone apparently due to the condensation 
of tetravalent counterions. At very high Cg (for instance, at Cg — 0.032 for (7+4 — 2.0), 
Q{r) in the zone can even attain a value smaller than the chain bare charge. In this case, 
negative particles are the major composition of the first ion layer around a chain, which is 
counterintuitive. The following three observations are worth to be noticed. (1) The salt 
concentration at which Q{r) becomes (roughly) zero in the region r > [1 + (T+4)/2 decreases 
with increasing (T+4. (2) The oscillatory behavior appears more markedly for large (74,4 than 
for small one. (3) At a given Cg, the larger the (7+4, the higher the maximum peak of Q{r) 
will be. 

It is tempting to consider that a polyelectrolyte and the ions condensed on the chain 
form a complex object. Under some conditions, the effective charge of the complex object 
reverses its sign. This phenomenon is called charge inversion [49] . Nguyen et al. intended to 
link the phenomenon of charge inversion with the reentrant condensation [5]. They proposed 
that zero effective charge and charge inversion of chains are, respectively, the cause for the 
condensation and for the redissolution of polyelectrolytes. Their theory necessitates further 
verification either by experiments or by simulations. In the rest of this section, we try 
to calculate the effective charge of a polyelectrolyte, in hope of being able to clarify the 
relationship between the charge inversion and the reentrant condensation. We simply define 
the complex to be the set of particles lying within a worm-like tube of radius r,. around 
the chain. This definition has been utilized in many places, for example, in Rcfs. [18] and 
[73] . Usually, Vc is chosen to be the distance at which the electrostatic attraction between a 
monovalent counterion and a monomer is equal to the kinetic energy ^k^T; in consequence. 
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the counterion inside the tube of radius rc can not escape from the chain. We, thus, consider 
the total charge Q{rc) inside the tube as the effective charge of a chain. In our study, Vc is 
equal to 2. The results are presented in Fig. 12(a). Two behaviors were observed: (1) for 
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FIG. 12: (a) Qir^) as a function of Cg with = 2.0, and (b) Q' as a function of Cg, for different 
C7+4. The symbol used to denote each c'"_|_4 is indicated in the figures. 



(T+4 = 3.0 and 4.0, Qijc) decreases with increasing Cg, (2) for cr+4 < 2.0, Qijc) increases 
basically with Cg. Notice that Qivc) for (74.4 = 1.4 or 2.0 attains a large positive value in 
the mid-salt region. Hence, IFmm'(^) should show strong repulsion in this region, according 
to the theory of Nguyen et al. But, in fact, we observed a strong attractive well for IFmm'('") 
and the chains were observed to aggregate together. Therefore, the results do not support 
the theory. 

We notice that it is too simple to regard Q{r-c) as the effective chain charge, because 
ion valence is not considered in the definition. For example, the suitable for tetravalent 
counterions should be 8.0, based upon the same argument to choice Tc for monovalent coun- 
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terions. In addition, Tc — 2.0 becomes too small while (7+4 > 3.0, since the depletion zone 
around a chain for the tetravalent counterions is bigger than the tube. Therefore, we rede- 
fine Tc = 8.0 as the novel condensation region for tetravalent counterions, whereas keeping 
Tc = 2.0 for monovalent counterions and for monomers on the other chains. Considering 
that the condensation of colons is mainly mediated by tetravalent counterions, we choose 
Tc — 8.0 for colons. The charges inside the different tubs are computed and their sum Q' 
is then regarded as an alternative definition for the effective chain charge. The results are 
presented in Fig. 12(b). We remark that Q' looks quite differently from Q{r^. It shows a 
positive peak in the mid-salt region for 17+4 = 1.0, 1.4, and 2.0 but retains negative value 
for the other (7+4. The obtained results suggest a disconnection between neutralization of 
chain charge and condensation of polyelectrolytes. 

Our simulating system in the high-salt region reveals a picture more close to the one 
described by Soils [52] : although locally overcompensated by condensed counterions, the net 
charge of a redissolved chain can be positive or negative due to the association of colons. 
We must point out that the effective chain charge interpreted by Q{rc) or by Q' in this 
section can be only served as a reference because both Q{r^ and Q' are sensitive to the 
choice of A rigorous way to study it, which can circumvent the ambiguity, is to study 
the electrophoretic mobility of polyelectrolytes in salt solutions. This part of research is 
currently undergoing. 

IV. CONCLUSIONS 

We have performed molecular dynamics simulations in canonical ensemble to investigate 
the properties of highly- charged polyelectrolytes in tetravalent salt solutions, including chain 
morphology, single-chain structure factor, swelling exponent, persistence length, like-charge 
attraction between chains and integrated charge distribution around a polyelectrolyte. We 
have explored a wide range of salt concentration and studied the effect of size of tetravalent 
counterions on the above properties. 

We have found that a polyelectrolyte exhibits, in turns, two conformational transitions 
upon addition of salt: a collapsed transition and a reexpanding transition, manifested by a 
V-shaped curve of the square of radius of gyration (i?g) against salt concentration Cg. The 
systems at different monomer concentration Cm show similar curves. The curves overlap 
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each other in the chain-reexpansion region. In the chain-collapsed region, we have also 
observed that AmdRl) ~ i^g^)) against Cg/Cm lie over each other for different Cm, where 
and (Rg^) weakly depends on Cm- These findings are useful in understanding the phase 
boundaries of the condensation of polyelectrolytes: on the upper boundary, Cg is roughly 
independent of Cm and on the lower one, Cg depends linearly on Cm- 

The study of chain morphology shows that for the cases with cr+4 (the size of the tetrava- 
lent counterions) around 1.0, the chains reveal compact and prolate structures in the vicinity 
of the equivalence concentration C*. In a low-salt region and in a high-salt region, the chains 
show extended structures. At a very high salt concentration, the chains behave similarly to 
neutral polymers. For the cases with large (T+4, the chains are elongated and do not display a 
coUapsed-transition in the mid-salt region. For the case with vanishing (T+4, the behavior of 
the chains is analogous to that of neutral polymers as Cg > C*. Noticeably, our results show 
that (7+4 = 1.0 is the optimal condition to pack polyelectrolytes into the smallest volume. 

We have shown that the single-chain structure factor S{q) of the multi-chain system in 
salt solutions follows the Guinier function A^(l — g^(i?g)/3) while q <IC R~^. S{q) then 
exhibits a power- law- like behavior q^ in the region R~^ <^ q <^ suggesting that the 
scaling law, R^ ~ N'^, is held for polyelectrolytes in salt solutions. The swelling exponent 
1/ — —1/s attains a minimum value near Cg* and behaves in analogy of the morphological 
quantities. It shows once more that the chains are mostly collapsed while (7+4 is compatible 
with 1.0. 

Our simulations incorporated explicitly the salt ions and, therefore, were able to study the 
negative regime of the electrostatic persistence length ig. This regime appears simultaneously 
with the mechanical instability of the polyelectrolytes induced by salt. We have examined 
the theory of Ariel and Andelman [48] and found that it fails to predict £e in the so- claimed 
valid region b <^ <^ {N — l)b but describes £e well in the unexpected region <^ b. 
Since the OSF theory is vahd in the latter region, we conjectured that the theory of Ariel 
and Andelman is applicable in the latter region. This part deserves a detail investigation in 
the future. 

We have studied like-charge attraction between polyelectrolytes by calculating the poten- 
tial of mean force VFmm' (f) between monomers on different chains. The results show that it 
depends strongly on salt concentration and ion size. An attractive region appears in Waim'ir) 
only when the salt concentration is intermediate and the tetravalent counterions possess a 
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size compatible with the monomers. The most profound attractive well is happened while 
(7+4 = 1.0 and Cs = C*. 

In the final part, the integrated charge distribution Q{r) around a polyelectroljd;e was 
investigated. We have found that the chain charge can be locally overcompensatcd by 
condensed ions. Q{r) displays an oscillatory behavior at high salt concentrations while 
c+4 > 1.0. The oscillatory behavior results from a multi-layer organization of ions around 
the chain, which is in agreement with the picture predicted by Solis [52]. We have used 
two different definitions to calculate the effective chain charge. The results suggest that 
charge inversion is not a necessary condition to occur reentrant condensation. Recently 
Wen and Tang have provided evidence against the connection between charge inversion and 
resolubilization of polyelectrolyte bundles [74]. A more rigorous study by electrophoresis 
performed by molecular dynamics simulations is undergoing. 
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